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Abstract 

In  this  study,  large  eddy  simulation(LES)  of  transonic  flow  around  the  NACA  0012  airfoil  is 
performed  accounting  for  the  Leonard  stress  terms,  the  cross-stress  terms  and  the  subgrid- 
scalc(SGS)  Reynolds  stress  terms  as  the  scale-similarity  model  at  a free  stream  Mach  num- 
ber of  0.8,  a Reynolds  number  of  9x10°  and  an  angle  of  attack  of  2.26°.  An  upwind  finite 
volume  formulation  is  used  for  the  discretization  of  compressible  spatial-filtered  Navier- 
Stokes  equations.  To  exclude  excessive  numerical  damping  due  to  the  shock-capturing 
scheme,  a hybrid  method  which  uses  linear  combination  of  the  third-order  upwind  scheme 
and  the  TVD  scheme  is  employed.  To  reduce  the  total  number  of  grid  points,  zonal  em- 
bedded mesh  is  employed  in  the  present  LES  analysis,  in  which  a computational  domain 
is  decomposed  near  the  wall-boundary.  In  the  case  represented  here,  it  is  shown  that  the 
statistical  values  in  the  turbulent  boundary  layer  with  shock/turbulence  interaction  is  able 
to  be  estimated,  and  characteristics  are  clarified  on  the  statistic  of  the  turbulence. 

1.  Introduction 

The  analyses  of  unsteady  flow  on  aircraft  and  fluid  machinery  in  the  transonic  speed  region 
will  be  required  in  future.  With  the  recent  development  of  the  super-computer,  LES  will 
be  used  for  such  complicated  flow  field  in  future.  Since  the  present  LES  analysis  needs  a 
shock-capturing  scheme  in  the  transonic  region,  it  is  required  to  investigate  the  effect  of  the 
inherent  numerical  dissipation  on  LES.  Referring  to  third-order  accurate  shock-capturing 
scheme,  the  third  order  numerical  viscosity  is  dominant  in  the  smooth  region,  whereas 
the  first-order  viscosity  is  dominant  in  the  region  where  the  discontinuity  appears.  The 
numerical  dissipation  due  to  the  third-order  numerical  viscosity  is  investigated  by  applying 
the  fourth-order  central  scheme  with  the  third-order  numerical  viscosity  to  the  turbulent 
channel  flow'.  As  the  results,  the  third-order  viscosity  with  appropriate  coefficient  is  found 
to  be  substituted  in  the  Smagorinsky  model  on  its  SGS  dissipation  rate.  Then,  in  this  study, 
eddy  viscosity  SGS  model  is  not  used  explicitly,  and  the  third-order  accurate  scheme  is 
used  to  substitute  the  SGS  dissipation  w'ith  the  appropriate  artificial  viscosity  coefficient. 
The  analysis  is  carried  out  taking  the  sum  of  the  Leonard  stress  terms,  the  cross-stress 
terms,  and  the  SGS  Reynolds  stress  terms  as  the  scale-similarity  model  into  account  on 
the  curvilinear  coordinates.  Still,  the  inviscid  numerical  flux  is  estimated  by  combining 
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the  third-order  scheme  with  the  TVD  scheme  in  order  to  stabilize  the  solution  near  the 
discontinuities. 

A single  structured  grid  is  not  suitable  for  LES  at  flight  Reynolds  numbers,  since  the 
viscous  sublayer  needs  the  very  fine  resolution  in  all  directions(Chapman,  1979);  As  the 
result,  the  structured  grid  has  numerous  grid  points  throughout  the  domain.  On  the  other 
hand,  the  unstructured  grid  is  efficient  to  avoid  the  unnecessarily  grid  resolution  for  much 
of  the  computational  domain( Jansen,  1999).  The  reduction  of  total  number  of  grid  points 
can  be  realized  by  the  embedded  mesh  system.  In  this  study,  we  have  chosen  the  embedded 
mesh  system  to  solve,  because  of  the  following  points:  1)  the  number  of  grid  points  can 
be  reduced  in  the  spanwise  and  streamwise  direction,  2)  the  computational  efficiency  is 
almost  the  same  per  unit  cell  in  comparison  with  the  efficiency  in  a single  structured  mesh 
system,  and  3)  it  is  comparatively  effortless  to  improve  the  spatial  accuracy. 

2.  LES  Methodology 

In  this  section,  we  summarize  the  LES  methodology,  including  governing  equations,  SGS 
model,  descritizing  formula  and  numerical  implementation. 


2.1.  GOVERNING  EQUATIONS 


The  governing  equations  are  the  spatially  filtered  compressible  Navier-Stokes  equations. 
The  spatial  filtering  removes  the  small-scale  components  of  the  fluid  motion,  retaining  the 
unsteadiness  associated  with  the  large-scale  turbulent  motion.  For  an  arbitrary  function 
/,  / represents  spacc-filtcrcd  variable.  For  compressible  flows,  it  is  expedient  to  define 
the  Favre  filtered  variable  /.  Applying  these  definitions  to  the  compressible  Navier-Stokes 
equations,  we  obtain  the  filtered  governing  equations 


dp  dpui  _ 

dt  dxi 


(1) 


dpuj  + d(pv.jUj  + 6jjp  - 9ij) drjj 

dt  dxj  Xj 


de  dje+pjUj  -a,^  - + dj) 

dt  + dxj  xj 


P 


(4) 


where  Xi  represents  the  Cartesian  coordinates(i=l,  2,  3),  p is  the  mean  density,  tq  are  the 
Cartesian  components  of  the  filtered  velocity,  p is  the  mean  pressure,  Ty  is  the  subgrid 
scale  stress  tensor,  <7y  is  the  molecular  stress  tensor,  e is  the  filtered  total  energy  per  unit 
volume,  qj  is  the  subgrid  heat  flux,  and  dj  is  the  energy  diffusion  by  the  subgrid  stresses. 
We  neglect  the  summation  term  th  in  the  equation  of  state,  which  would  be  very  small 
compared  to  the  thermodynamic  pressure  p. 

2.2.  DESCRITIZATION 

The  governing  equations  can  be  descritized  in  finite  volume  method  for  a control  volume 
V with  surface  dV : 

dt  Jv  4-  E2ri]2  + Ean\s)dA  = 0, 


(5) 


Best  Available  Copy 


LARGE  EDDY  SIMULATION  OF  TRANSONIC  TURBULENT  FLOW  745 


where 


p 

PUi 

<?  = 

pUi 

, Ei  = 

pUiiij  + Syp  - dij  4-  Ty 

e 

_ puiht  - - dijUi  + <j;  + d{ 

(ftiii fti20U3)#  denotes  the  normal  vector  to  each  cell-interface.  The  explicit  Runge  Kutta 
method  is  employed  for  the  time  integration.  The  third-order  numerical  inviscid  flux 
is  estimated  by  using  the  fourth-order  central  scheme  and  the  additional  third-order 
dissipation. 


r(.W)  _ 

•'i+i/2  ~ 

fi +1/2  + C'l(-ftl^!)i+l/2(Q!i+3/2  _ 2ai+i/2  + 1/2)1 

(7) 

11 

3 + 

1/2'  Qft  l/2)’ 

(8) 

Qi+ 1/2  — 

^i+l/2^i  + 1/2^’ 

(9) 

where  R and  A represent  the  diagonalization  matrix  of  the  right  eigenvectors  and  the 
diagonal  matrix  of  the  eigenvalues,  respectively.  The  values  of  each  variables  on  either  side 
of  cell  interfaces  are 

Q(+1/'2  = Qi+  fiQ  + 6Ai+1//2C?  - Ai+3/2Q), 

Qiii/2  = Qi  ~ + 6Ai_,/2Q  - Aj_3/2Q).  (10) 

The  values  of  artificial  viscosity  coefficient  £4  is  calibrated  at  1 /256  in  the  airfoil  simulation. 

The  method  mentioned  above  cannot,  treat  the  transonic  flow,  since  the  shock  wave 
appears.  Then,  the  numerical  flux  fj+yi 's  differently  prepared  by  using  a finite  difference 

split.ting(Roe,  1981)  and  extrapolated  variables  Q?±y2  with  a slope-liiniter(Daiguji  et 
al,  1997).  To  exclude  excessive  numerical  damping  due  to  the  limiter,  we  employ  self- 
adjusting  hybrid,  that  uses  linear  conbination  of  the  third-order  upwind  scheme  and  the 
TVD  scheme,  the  numerical  flux  /;+1/2  is  blended  as  follows: 


/i+j/2  — (1  ~ f?)/j+i/2^3r<^  + 0/i+ 1/2> 


(11) 


f 0.  for  < 0.01 

Q =:  ) p,+l+2p;-pi-l  “ 

1 1.  for  > 0.01 

^ ’ Pi+i +2p,4-pi_i 


(12) 


2.3.  SGS  MODEL 

SGS  stresses  can  be  decomposed  into  the  following  terms(Leonard,  1974): 

Ty  = Lij  + Cij  + Rij,  (13) 


the  Leonard  stress  terms,  Ly  = p(wjfy  — Ujfy);  the  cross-stress  terms,  Cij  — p(h,-u' 
and  the  SGS  Reynolds  stress  terms,  Rij  = pu'w'.  Present  study  follows  the  hypothesis  that 
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the  SGS  energy  dissipation  would  be  substituted  in  the  numerical  dissipation  due  to  the 
inherent  artificial  viscosity  of  the  shock-capturing  scheme.  Consequently,  the  eddy  viscosity 
model  was  not  employed  in  the  present  LES  analyses.  However,  a scale-similarity  model  is 
taken  into  account.  They  can  be  estimated  using  Taylor  expansion.  To  eliminate  the  terms 
which  violate  the  Galilean  invariance(Speziale,  1985),  these  terms  are  summated  with  the 
Bardina  model(Bardina  et  al,  1980)  being  applied  to  the  cross-stress  terms  and  the  SGS 
Reynolds  stress  terms. 


Lij  + Cij  + Rij  « p{uiUj-UiUj) 


(14) 


For  the  SGS  terms  qj  and  d3  in  energy  equation,  the  following  models  are  employed 
by  approximating  u3  = Uj  and  T = T. 


q3  = Cpp(ujT  — UjT)  = Cpp(ujT  - UjT) 
dj  - ~p(u2Uj  - u2Uj) 

+ Aj), 


24  k ~dxk  dxk  ' 2i“kdxkdxk 
where  a perfect  gas  with  a constant  specific  heat  capacity  Cp  is  assumed. 


(15) 


(16) 


3.  Results  and  Discussion 
3.1.  CODE  VALIDATION 


A validation  test  for  the  shock-capturing  is  performed  for  the  shock  tube  problem(Shu  et 
al. , 1989),  which  has  several  extrema  in  the  smooth  region.  Figure  1 shows  the  solution 
via  the  present  method  with  400  cells.  For  the  comparison,  we  have  included  the  line  of 
solution  via  the  8th-order  ENO  scheme(Harten,  1989)  with  1600  cells.  Note  that  there  is 
no  serious  overshooting  or  undershooting  neat  the  discontinuity.  This  fact  shows  that  the 
present  method  is  suitable  for  LES  with  strong  discontinuities. 

A validation  test  for  large  eddy  simulation  is  carried  out  by  applying  the  present 
method  to  the  channel  flow.  The  analyses  are  performed  for  the  turbulent  plane  Couette 
flow  of  8.0/ix2.0hx4.0/i  at  Rc(3 * * &^)  = 1300  and  the  Mach  number  M(^)  = 0.3.  The 
number  of  cells  is  (64,64,32)  with  the  size  of  cells  in  the  wall  coordinate  varying  from  1 
to  1.4  wall  units.  Figure  2 shows  the  mean  velocity  profiles  via  the  TVD  scheme  and  the 
present  method.  The  ensemble  data  are  obtained  from  the  averaged  flow  over  100000  time 
steps.  We  have  included  incompressible  DNS  results(Bech  et  al,  1995)  for  the  comparison. 
The  artificial  coefficient  £4  is  calibrated  at  1/256,  that  is  also  used  in  the  present  airfoil 
simulation.  In  this  case,  the  hybridization  of  the  scheme  is  absolutely  necessary  since 
the  inherent  numerical  dissipation  due  to  the  flux-limiter  of  the  TVD-scheme  becomes 
dominant. 

In  Fig.  3 we  show  the  mean  velocity  profiles  with  and  without  the  SGS  models  for 
the  same  channel  flow.  We  have  included  the  root  mean  square  ( mis ) of  the  velocity 
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components  in  Fig.  4.  The  predictions  in  all  rms  quantities,  especially  the  streamwise 
component,  are  of  great  improvement  by  using  the  present  SGS  models,  although  a little 
hit  greater  artificial  coefficient.  S4  of  1/128  is  taken  to  produce  more  dissipative  solution 
compared  to  the  DNS  results. 

3.2.  AIRFOIL  SIMULATION 

The  flow  configuration  is  that  of  experiments(Harris,  1981)  at  a Reynolds  number  based 
on  chord  lie  = 9 x 106,  a free -stream  Mach  number  M — 0.8,  and  an  angle  of  attack 
o = 2.2GD . A single  structured  grid  system  makes  it  difficult  to  perform  LES  for  the  flow 
of  the  present  Reynolds  number.  The  use  of  embedded  mesh  reduces  the  computational 
memory  and  the  CPU  to  some  extent.  The  mesh  size  is  refined  near  the  wall  and  the  mesh 
system  is  replaced  6 times  toward  the  wall  surface  regarding  the  streamwise  direction 
and  the  spamvise  direction  so  that  the  large  eddy  can  be  resolved,  as  shown  in  Fig.  5. 
The  extent,  of  spanwise  direction  is  chosen  as  Lz  — 0.04c.  Current  simulation  employs 
approximately  11  million  cells. 

Figure  6 shows  the  instantaneous  Mach  number  distribution  in  the  mid-span  plane. 
On  the  suction  side,  the  shock  wave  appears  at  about.  60%  of  the  chord  length  from  the 
leading  edge.  It  is  found  that  there  is  no  serious  numerical  oscillation  near  the  shock  wave. 
Figure  7 shows  the  isosurfaces  of  instantaneous  spanwise  velocity  w on  the  suction  side. 
Weak  three  dimensionalization  occurs,  and  secondary  flow  is  formed  (Fig.  8).  It  is  also 
shown  that  the  three-dimensionalized  structure  becomes  larger,  as  it  goes  to  the  trailing 
edge,  as  shown  in  Fig.  9.  The  increase  of  the  friction  coefficient  can  be  also  confirmed  with 
the  relation  to  the  three-dimensionalization.  In  Fig.  10,  we  show  the  friction  cocficient, 
which  increases  at  about  5%  of  the  chord  length  from  the  leading  edge  on  the  suction  side 
where  the  transition  takes  place  naturally. 

Figure  11  show's  mean  wall  pressure  distribution,  which  almost  agrees  with  the  ex- 
perimental result,  however  there  is  a little  discrepancy  of  the  shock  position  between  the 
simulation  and  the  experimental  data.  The  aspect  of  the  pressure  distribution  under  the 
influence  of  adverse  pressure  gradient  by  the  shock  -wave  is  predicted  so  as  to  produce  the 
B.  L.  separation  more  extensively. 

Figures  12  and  Fig.  13  show  the  turbulence  energy  profiles  on  the  pressure  side  and  the 
suction  side,  respectively.  The  djst.ance  from  the  wall  is  represented  by  rj.  On  the  suction 
side,  the  turbulence  energy  rises  (x/c  — 0.7)  behind  the  shock  wave,  while  it  contrastively 
decreases  downstream  the  position  of  (x/c  — 0.4)  on  the  pressure  side. 

The  spamvise  one  -dimensional  energy  spectrum  distribution  is  examined  in  order  to 
understand  whether  it.  reproduces  the  energy  cascade  process  in  the  turbulent  boundary 
layer.  The  energy  spectrum  is  defined  as, 

Eiti(k)  = - / Ri^e-^dr,  (17) 

7T  JO 

where  Rhi  represents  the  correlation  tensor  given  by 

Ri,i(r)  =<  u'iWu'ii x + r)  > . (18) 

v(  is  approximated  by  using  the  grid-scale  differences  u,—  < u,  >.  r is  in  spanwise  direction 
in  this  case.  Figure  14  shows  the  spanwise  one-dimensional  energy  spectra  of  the  velocity 
at  y+  ~ 10.  It  is  confirmed  that  the  energy  cascade  with  the  close  gradient  of  the  5/3  law 
in  the  sampling  point  (x/r.  = 0.4, 0.7)  is  obtained. 
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4.  Conclusions 

The  scale-similarity  models  are  taken  into  account  in  the  generalized  curvilinear  coordinate 
system  in  order  to  represent  the  SGS  effect  which  includes  the  backscatter  in  the  energy 
cascade  in  turbulent  boundary  layer.  To  keep  the  resolution  of  the  near  wall  turbulence, 
the  embedded  mesh  system  is  employed  in  this  study.  Also,  an  highly  accurate  scheme  is 
obtained  in  order  to  capture  generation  of  the  turbulent  boundary  layer  from  the  leading 
edge  and  its  developing  process.  Concurrently,  the  present  hybridization  of  the  scheme 
produces  no  serious  oscillation  near  the  shock  wave  in  the  present  simulation.  As  a result 
of  applying  the  present  method  to  transonic  turbulent  flow  around  the  airfoil,  it  is  possible 
to  reproduce  the  flow  field  with  large-scale  structure  of  the  turbulent  boundary  layer.  It 
is  found  that  the  qualitative  analysis  is, possible  on  the  turbulent  characteristics,  that 
includes  the  interaction  with  shock  wave  by  the  present  LES  method. 
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Figuiv  1.  Code  verification  for  the  shock-  cap- 
turing 


Figure  5.  Side  view  of  embedded  mesh  around 
the  NACA0012  airfoil 


Figure  2.  Code  verification  for  the  channel  tur- 
bulent (low 


Figure  3.  Model  verification  for  the  channel 
turbulent  flow 


Figure  Model  verification  for  the  channel 
turbulent  flow 


Figure  6.  Side  view  of  instantaneous  Mach 
number  contours 


Figure  7.  Isosurfaces  of  spanwise  velocity 


Figure  8.  Isosurfaces  of  spanwise  velocity  near 
the  leading  edge 
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Figure  9.  Isosurfaccs  of  spanwisc  velocity  near 
the  trailing  edge 


Figure  12.  Turbulence  energy  profiles  on  the 
suction  side 


Figure  10.  Skin  friction  distribution  Figure  13.  Turbulence  energy  profiles  on  the 

pressure  side 
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Figure  11.  Wall  pressure  distribution 


Figure  14-  Spanwise  one  dimensional  energy 
spectra  on  the  pressure  side  at  y+  2;  10 
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